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Abstract 

We investigate a model for traffic flow based on the LighthiU-Whitham-Richards model 
that consists of a hyperbolic conservation law with a discontinuous, piecewise-linear flux. 
A mollifier is used to smooth out the discontinuity in the flux function over a small 
distance e ^ 1 and then the analytical solution to the corresponding Riemann problem 
is derived in the limit as e — > 0. For certain initial data, the Riemann problem can 
give rise to zero waves that propagate with infinite speed but have zero strength. We 
propose a Godunov-type numerical scheme that avoids the otherwise severely restrictive 
CFL constraint that would arise from waves with infinite speed by exchanging informa- 
tion between local Riemann problems and thereby incorporating the effects of zero waves 
directly into the Riemann solver. Numerical simulations are provided to illustrate the be- 
haviour of zero waves and their impact on the solution. The effectiveness of our approach 
is demonstrated through a careful convergence study and comparisons to computations 
using a third-order WENO scheme. 
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1. Introduction 

In the 1950's, Lighthill and Whitham 30J and Richards JB] independently proposed 
the first macroscopic traffic flow model, now commonly known as the LWR model. Al- 
though this model has proven successful in capturing some aspects of traffic behaviour, 
its limitations are well-documented and many more sophisticated models have been pro- 
posed to capture the complex dynamics and patterns observed in actual vehicular traffic 
[22] . Despite this progress, the LWR model remains an important and widely- used model 
because of its combination of simplicity and explanatory power. 
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The LWR model consists of a single scalar nonlinear conservation law in one dimension 



Pt + f{p). - 0, 



(1) 



where p{x,t) is the traffic density (cars/m), 

fip) = pv{p) 

is the traffic flow rate or flux (cars/sec), and v{p) is the local velocity (m/sec). The most 
commonly used flux function is 



which was obtained by Greenshields in the 1930's by fitting experimental measure- 
ments of vehicle velocity and traffic density. Here, Umax is the maximum free-flow speed, 
while Pmax is the maximum density corresponding to bumper-to-bumper traffic where 
speed drops to zero. The LWR model belongs to a more general class of kinematic wave 
traffic models that couple the conservation law Eq. ([I]) with a variety of different flux 
functions. 

Extensive studies of the empirical correlation between flow rate and density have 
been performed in the traffic flow literature. This correlation is commonly referred to 
as the fundamental diagram and is represented graphically by a plot of flux / versus 
density p such as that shown in Fig. [T] A striking feature of many experimental results 
is the presence of an apparent discontinuity that separates the free flow (low density) 
and congested (high density) states, something that has been discussed by many authors, 
including [3 [TH [TSl [53]. In particular, Koshi et al. [57 characterize flux data such as 
that shown in Fig. [T] as having a reverse lambda shape in which the discontinuity appears 
at some peak value of the flux. 

This behavior is also referred to as the two-capacity or dual-mode phenomenon [31 13] 
and has led to the development of a diverse range of mathematical models. Zhang and 
Kim [41] incorporated the capacity drop into a microscopic car-following model that 
generates fundamental diagrams with the characteristic reverse-lambda shape. Wong 
and Wong ,38J performed simulations using a multi-class LWR model from which they 
also observed a discontinuous flux-density relationship. Colombo 8 and Goatin [17] 
developed a macroscopic model that couples an LWR equation for density in the free 
flow state, along with a 2x2 system of conservation laws for density and momentum in 
the congested state; the phase transition between these two states is a free boundary 
that is governed by the Rankine-Hugoniot conditions. Lu et al. 1311 incorporated a 
discontinuous (piecewise quadratic) flux directly into an LWR model, and then solved 
the corresponding Riemann problem analytically by constructing the convex hull for a 
regularized continuous flux function that consists of two quadratic pieces joined over a 
narrow region by a linear connecting piece. 

There remains some disagreement in the literature regarding the existence of dis- 
continuities in the traffic flux, with some researchers (e.g.. Hall \2T) arguing that the 
apparent gaps are due simply to missing data and can be accounted for by providing 
additional information about traffic behaviour at specific locations. Indeed, Persaud and 
Hall [3? and Wu [3^ contend that the discontinuous fundamental diagram should be 
viewed instead as the 2D projection of a higher dimensional smooth surface. 




(2) 
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We will nonetheless make the assumption in this paper that the fundamental diagram 
is discontinuous. Our aim here is not to argue the validity of this assumption in the 
context of traffic flow, since that point has already been discussed extensively by [311 
[551 EI], among others. Instead our objective is to study the effect that such a flux 
discontinuity has on the analytical solution of a ID hyperbolic conservation law, as well 
as to develop an accurate and efficient numerical algorithm to simulate such problems. 

A related class of conservation laws, in which the flux /(p, x) is a discontinuous func- 
tion of the spatial variable x, has been thoroughly studied in recent years (see [SlIS] and 
the references therein) . Considerably less attention has been paid to the situation where 
the flux function has a discontinuity in p. Gimse [16| solved the Riemann problem for a 
piecewise linear flux function with a single jump discontinuity in p by generalizing the 
method of convex hull construction [551 Ch. 16]. In particular, Gimse identified the exis- 
tence of zero shocks, which are discontinuities in the solution that carry no information 
and have infinite speed of propagation. We note that more recently, Armbruster et al. [T] 
observed zero rarefaction waves with infinite speed of propagation in their study of supply 
chains with finite buffers (although they did not refer to them using this terminology). 

Gimse's results were improved on by Dias and Figueira [11^. who used a mollifier 
function 77^ to smooth out discontinuities in the flux function over an interval of width 
< e <C 1 before constructing the convex hull using standard techniques. Solutions to 
the mollified problem were proven to converge to solutions of the original problem in the 
limit as e — > [11 . Dias and Figueira's framework has also been applied to problems 
involving fluid phase transitions [TOl [13] and viscoelasticity [12^ . 

In this paper, we apply Dias and Figueira's moUiflcation approach to solving a conser- 
vation law with a piecewise linear flux function f(p) in which there is a single discontinuity 
at p = pm (see Fig.jl]). The model equations and their relevance in the context of traffic 
flow are discussed in Section [2j We introduce a mollified flux function /e(p) in Section 
[3] and verify its convexity, which then permits us to derive the analytical solution to the 
Riemann problem using the method of convex hull construction. 

In Section |4| we consider the special case where either of the two constant initial 
states in the Riemann problem equals p,„, the density at the discontinuity point. This 
is precisely the case when a rarefaction wave of strength 0(e) and speed 0(l/e) arises, 
which approaches a zero rarefaction in the limit of vanishing e. There are two issues that 
need to be addressed regarding these zero waves. First, we consider the convergence of the 
mollified solution to that of the original problem, since Dias and Figueira's convergence 
results llj do not consider (nor easily extend to) the case when the left or right initial 
states in the Riemann problem are identical to pm- Secondly, we discuss the physical 
relevance of an infinite speed of propagation in the context of traffic fiow. 

The remainder of the paper is focused on constructing a Riemann solver that forms 
the basis for a high resolution finite volume scheme of Godunov type. Because zero waves 
travel at infinite speed, the usual CFL restriction suggests that choosing a stable time step 
might not be possible. Some authors have avoided this difficulty by using an implicit time 
discretization |33| . but this approach introduces added expense and complication in the 
numerical algorithm. Another approach employed in [l] is to replace the discontinuous 
flux by a regularized (continuous) function which joins the discontinuous pieces by a linear 
connection over an interval of width e ^ 1, after which standard numerical schemes may 
be applied; however, this approach requires a small e to achieve reasonable accuracy 
resulting in a severe time step restriction. 
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Figure 1: In |(a)| the discontinuous "reverse lambda" flux function Eq. (j4j is overlaid with empirical 
data extracted from 1211 Fig. 7] (reproduced with permission of Elsevier B.V.), along with the quadratic 
Greenshields flux The mollified flux from Eq. ^ is depicted in |(b)[ 



94 We use an alternate approach that ehminates the severe CFL constraint by incorpo- 

95 rating the effect of zero waves directly into the local Riemann solver. In the process, we 

96 find it necessary to construct solutions to a subsidiary problem that we refer to as the 

97 double Riemann problem, which introduces an additional intermediate state correspond- 

98 ing to the discontinuity value p = pm- A similar approach was used by Gimse |16j who 

99 constructed a first-order variant of Godunov's method, although he omitted to perform 

100 any computations using his proposed method. We improve upon Gimse's work in three 

101 ways: first, we solve the double Riemann problem within Dias and Figueira's moUifica- 

102 tion framework; second, we implement a high resolution variant of Godunov's scheme to 

103 increase the spatial accuracy; and third, we provide extensive numerical computations 

104 and a careful convergence study to demonstrate the effectiveness of our approach. 



2. Mathematical Model 

We are concerned in this paper with the scalar conservation law 

Pt + fip). = 0, 

having a discontinuous flux function 



fip) 



gf{p), if ^ p < prn, 

9c{p), ifpm^ps^l, 



(3) 



(4) 



that is depicted in Fig. l|a) The vehicle density p{x,t) is normalized so that ^ p ^ 1, 
and Pm is the point of discontinuity in the flux f{p). We restrict the flux to be a piecewise 



linear function in which the free flow branch has 
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9f{p) = P, 

and the congested flow branch has 

9c{p) = 7(1 - P)- 

10s Experimental data suggests that gf{pm) > gciPm), and so we impose the constraint 

109 We utilize the moUifier approach of Dias and Figueira 1,1 11 in order to approximate 
no the original equation ([S]) by 

111 where the mollified flux (pictured in Fig. |l|[b)[ ) is 

MPe) = Pe + il - h + ^)Pe) [ T],{s - p^n) ds, (7) 

112 with < e <C 1. The moUifier function is given by rj^{s) = ^ri{s/e), where ry(s) is a 

113 canonical moUifier that satisfies the following conditions: 

(i) ?7 ^ 0; 

115 (ii) rj E C°°(M) and is compactly supported on [—1, 1]; 

116 (iii) s) — ri{s) for all s G M; and 

The results in Dias and Figueira Jill guarantee that any moUifier satisfying the above 
criteria converges to the same unique solution in the limit e — )■ 0. We use the following 
moUifier 

f Cexp(^), if|s|<l, 
1 0, if |,s| ^1, 

121 where C w 2.2522836 ... is a constant determined numerically so that condition ( pv| ) 

122 holds; this choice is made for reasons of analytical convenience since the derivative rj' can 

123 be written in terms of rj. 
Because the mollified fiux function is smooth, the conservation law ^ may now 

be solved using standard techniques. We note that in the context of traffic flow, a 
potential problem arises when applying the usual Olemik entropy condition |34| as the 
selection principle to enforce uniqueness. Although Olemik's entropy condition does 
yield the physically-correct weak solution in the context of fluid flow applications, it 
does not always do so for kinematic wave models of traffic flow (see LeVeque [37], for 
example) . In particular, applying Olemik's entropy condition can lead to a solution that 
is not anisotropic [9] , corresponding to the non-physical situation where drivers react to 
vehicles both in front and behind. 
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Zhang PO' suggests two additional criteria on the flux function to guarantee anisotropic 
flows in kinematic wave traffic models. First, the characteristic velocity should be smaller 
than the vehicle speed v{p) = f{p)/p. That is, we require 

^>f{p), VO<p<l. 
P 

Secondly, all elementary waves must travel more slowly than the vehicles carrying them, 
or 

fiPl) f{Pr)\ ^ I{Pr)-f{pi) 



Pi Pr J Pr - Pi 

for all < Pr < 1- As shown in [57], our flux function Q satisfies both of these 
conditions and therefore it is reasonable to apply the Olemik entropy condition as the 
selection criterion for our traffic flow model. 



3. Exact Solution of the Riemann Problem with MolUfied Flux 
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We next construct and analyze the solution Pe(a;, t) of the mollified Riemann problem, 
which consists of the conservation law ^ along with flux ^ and piecewise constant 
initial conditions 

pi, ifx<0, 

Pr, if X ^ 0. 

This problem can be solved using the method of convex hull construction [551 Ch. 16] 
which requires knowledge of the inflection points of the mollified flux f^(p^). Using 
Eq. ([8| , the first and second derivatives of the flux are 



Peix,0) 



Ki.P<i) = 1 + [7 - (7 + VtiPe - Pm) - (7 + 1) / V^is - Pm) ds, 

2" 



fJip,) = K{p,) 



7 



7 + 1 



Pe [Pm ~ Pe) 



Pe - Pn 



1 



(9) 

(10) 



where 



2?7e(p<: - Pm) 



Since the moUifier has compact support, we know that there can be no infiection points 
outside the smoothing region of width 2e; that is, f'J{pt) = when \p^ — Pm\ ^ £■ Also, 
since K{pe) > 0, the convexity of is determined solely by the sign of the quantity 

-I 2 



PiPe 



which we analyze next. 
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7+1 



(Pn 



Pn 



1 



(11) 



Since Eq. ( 11 ) is a quartic polynomial in p^, analytic expressions are available for the 



roots; however, these are too complicated for our purposes. Instead, we take advantage 
of the scaling properties of the polynomial to simplify P and rewrite Eq. ( 11 ) as 



' 6 



(12) 



144 where y = Pe - Pm & [-e, e] and M = p,„ -7/(7 + 1). Clearly \y\ = \pe - Pm\ < e, and 

145 SO it is natural to suppose that y = 0{£^) with A > 1, after which the various terms in 



146 Eq. (12 1 have the scalings indicated below: 

1 
e 

0(e4A-2) o(e2A) O(e^) 0(£2) 



piy) - -72^' + '^y' + - ■ (13) 



Since M is a constant that is independent of e, we can take M = 0(1) as e — ?> 0. 



Guided by the scalings in Eq. (13 1, the dominant terms in P{y) are the last two terms 
having orders 0{e^) and O(e^). When e is sufficiently small, we may therefore neglect 
the remaining terms and determine the convexity of based on the sign of the simpler 
linear polynomial 

Piiy)^My-e^, 

which has a single root at y = /M corresponding to pe = Pm + e^/M. For values of 
y close enough to e^/Af , higher order terms in the polynomial P{y) become significant 
and could potentially introduce additional roots; however, by continuing this method of 
dominant balance, we find that P{y) maintains the single root at higher orders as well, 
which we demonstrate numerically using the plots of P summarized in Fig. [2j Based on 
this argument and the observation that inequality ([5| requires M > 0, we can conclude 
that 

/"(pe) ^ when p^ £ [p„,,, 1], 
and /"(pe) ^ when peG[0,Pm]. 

Therefore, the mollified flux function has a single inflection point at Pe — p„j + O(e^) as 
e — 0, where the slope /^(pm) — > —00. Since the flux derivative /'(p) corresponds to the 
elementary wave speed in the Riemann problem, this same point pm is also the source 
of the infinite speed of propagation which will be analyzed in more detail in Section [4] 

Using this information, we can now construct the convex hull of the flux function 
fe{Pe) which is then used to solve the Riemann problem. There are three non-trivial 
cases to consider, depending on the left and right initial states, pi and pr- Two of these 
cases (which we call A and B) lead to the emergence of a new constant intermediate 
state with density pm- This "plateau" is a characteristic feature of solutions to our LWR 
model with discontinuous flux. 

Case A: Pr < Pm < Pi- 

Here we construct the smallest convex hull of the set {(p^,y) : pr < Pe < pi and y ^ 



161 



/c(pe)}, which as shown in Fig. i 'a) must consist of three pieces. The first piece corre- 
sponds to a contact line that follows /e(pe) on the left up to the point p* G {pm — e, Pm) 
for which the shock and characteristic speeds are equal; that is, 

MPl)-fe{p*) . 
Pi - P* 

The third piece of the convex hull corresponds to a shock that connects the states p* and 
p/. The middle piece, in between the contact line and shock, gives rise to a rarefaction 
wave that follows the curved portion of the flux in the neighbourhood of pm- Based on 
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Figure 2: Plots of the function P(ez) and its roots for different e when 7 = 0.5 and = 0.5. 



this convex hull, we can then construct the solution profile shown in Fig. |t|[b)[ Since the 
rarefaction wave consists of density values bounded between pm — e and pm + e, this wave 
flattens out and degenerates to a constant intermediate state p„ in the limit as e — 0. 
To summarize, in the limit as e — > the solution to the Riemann problem when 

Pr < Prn < Pi IS 

p;, if a; < st, 

p{x,t)^{ p,n, ifst^x^t, (14) 

Pr, if X > t, 



consisting of a 1-shock moving to the left with speed 

„_ f{pi)-9f{Pm) 



Pi - Pm 

and a 2-contact moving to the right with speed 1. 



< 0, (15) 



Case B: 7/ (7 + 1) < < p™, < Pr- 

In contrast with the previous case, we now have pi < pr and so the convex hull 
: pi < pe < Pr and y ^ /^(pe)} lies above the flux-density curve as shown in 
Fig. ^[a) The first piece of the convex hull corresponds to a shock wave connecting the 



states pi and G {pm — £,Pm), while the third piece follows the portion of the flux 
function in the region [p^,,pr]. As before, the state p^ is chosen so that the shock and 
characteristic speeds are equal: 

Mp*) - fejpi) f,, . 

s = = hip*)- 

p* - Pl 
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Figure 3: The convex hull construction for Case A. |(a)| The thin (blue) line is the mollified flux function 
and the thick (black) line is the convex hull. The shaded oval region highlights the rarefaction 
wave. 1(b) I The solution p^(x/t) is shown as a thick (black) line, along with the corresponding flux 
derivative /'(p) as a thin (blue) line. Parameter values are p; = 0.9, pr = 0.2, pm = 0.5, e = 0.1 and 
7 = 0.5. 



171 The resulting solution profile in the limit as e — > is pictured in Fig. |^|[^b)| 

{pi, ii X < st, 
Pm, if st < a; ^ -7i, (16) 
Pr, if a; > -7i, 

172 which consists of a 1-shock moving to the left with speed 

gc{Pm) - f{pi) 



173 
174 
175 
176 
177 
178 
179 



Pm - Pi 

and a 2-contact moving to the left with speed —7. 



< 0, (17) 



Case C: pi < Pm < pr and 7/(7 + 1) ^ p;. 

The solution structure in this case is significantly simpler than the previous two in that 
there is only a single shock connecting the states p^. and pi, and hence no intermediate 
state. The convex hull is depicted in Fig. I a) and the corresponding solution profile in 
Fig. |^[b)[ As e — > 0, the solution reduces to 

, s ( Pi, if a; < st, , , 

= \ Pr, ifx^st, (1^) 

which corresponds to a shock with speed 

^^ fiPr)-f{pi) ^ (19) 

Pr ~ Pi 
9 




Figure 4: The convex hull construction for Case B. |(a)| The thin (blue) line is the mollified flux function 
fe{Pi) and the thick (black) line is the convex hull. The shaded oval region highlights the rarefaction 
wave. 1(b) I The solution p^{x/t) is shown as a thick (black) line, along with the corresponding flux 
derivative /'(p) as a thin (blue) line. Parameter values are p; = 0.4, pr = 0.9, e = 0.1, 7 = 0.5, 
Pm = 0.5. 




Figure 5: The convex hull construction for Case C. |(a)| The thin (blue) line is the mollified flux function 
feiPa) and the thick (black) line is the convex hull. The shaded oval region highlights the rarefaction 
wave. 1(b) I The solution p^(x/t) is shown as a thick (black) line, along with the corresponding flux 
derivative /'(p) as a thin (blue) line. Parameter values are p; = 0.2, pr = 0.9, pm = 0.5, e = 0.1 and 
7 = 0.5. 
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Figure 6: Four possible cases for the initial conditions of the double Riemann problem |[20l 



which can be either positive or negative depending on the sign of the numerator. 

Note that as pi or approaches the discontinuity Pm, the Riemann solution becomes 
sensitive to the choice of initial data. This sensitivity is not unique to our problem, but is 
also observed in other analytical solutions for problems with discontinuous flux, such as 
in [13l [TBI [31] • In the context of traffic flow, this sensitivity occurs in the neighbourhood 
of the transition point between free-flow and congested traffic, the exact location of which 
is expected to be highly sensitive to the state of individual drivers comprising the flow. 
Therefore, the sensitivity in our model is consistent with actual traffic. 



189 4. Analysis of Zero Waves 

190 In this section, we consider the two special situations that were not addressed in 

191 Section [Sj namely where either pi = pm or pr — Pm- In both cases, the mollified problem 

192 gives rise to a wave having speed 0(l/e) and strength 0(e), which we refer to as a zero 

193 rarefaction wave because of its similarity to the zero shocks identified by Gimse |16j . 

194 Since the speed of these waves becomes infinite as e — 0, information can be exchanged 

195 instantaneously between neighbouring Riemann problems in any Godunov-type method. 

196 We demonstrate in this section how these effects can be incorporated into the local 

197 Riemann solver. 

198 Motivated by the need to consider interactions between Riemann problems arising 

199 from two pairs of piecewise constant states, we consider a double Riemann problem con- 

200 sisting of two "usual" Riemann problems: one on the left with (pi,pr) — {Ci,pm), and a 

201 second on the right with (pi,pr) = (prmCr)- As a result, the mollified conservation law 

202 ^ is supplemented with the following piecewise constant initial data pictured in Fig. [6] 

{C'l, ifa;<a;i, 
Pm, if xi 5$ a; < X2, (20) 
Cr, if X > X2, 

203 where we have used the notation Ci and Cr for the left /right states to emphasize the 

204 fact that we are solving a double Riemann problem. 

205 Since the solution can contain waves whose speed becomes unbounded, it follows 

206 that we cannot solve the double Riemann problem with intermediate state Pm by simply 

11 



207 splitting the solution into two local Riemann problems and applying standard techniques. 

20S Instead, the origin and dynamics of zero waves need to be considered when constructing 

209 the solution. The double Riemann problem with a mollified flux consists of up to three 

210 separate waves. Two of the waves - which we call the 1- and 2- wave using the standard 

211 terminology - correspond respectively to the left-most (xi) and right-most (X2) waves 

212 arising from the initial discontinuities in the double Riemann problem. We will see later 

213 that the 1-wave is a shock and the 2-wave is a contact line. The third wave corresponds 

214 to a left-moving rarefaction wave having strength 0{e) and speed 0(l/e) that originates 

215 from the X2 interface and is located between the other two waves. As e — >■ 0, the 

216 rarefaction wave approaches a zero wave that mediates an instantaneous interaction 

217 between the 1- and 2-waves. 

21s 4.1. Origin of Zero Waves 

219 Since the local Riemann problem arising at the 2-wave (i.e., the contact line at the X2 

220 interface) is the source of the zero rarefaction wave, we begin by focusing our attention 

221 on the right half of the double Riemann problem. The formation of a zero wave at the 

222 discontinuity in the initial data located at X2 can be divided into two cases, corresponding 

223 to whether C'r < Pm or Cr > Pm- The specifics of the interaction between the 1-shock 

224 and the zero wave will be treated separately Section [4?2| 

225 

226 Zero rarefaction with Cr < Pm and Ci — Pm- 

227 We first consider the mollified Riemann problem with Ci = pi = pm and Cr = Pr < 

228 Pm- The convex hull and solution profile shown in Fig. [7] exhibit a right-moving contact 

229 line (the 2-wave) having speed s — 1 and a rarefaction wave of strength 0(e) travelling to 

230 the left with speed 0(l/e). For a simple isolated Riemann problem, the solution would 

231 reduce to a lone contact line as e — > and the zero rarefaction would have no impact. 

232 However, when the zero wave is allowed to interact with the solution of another neigh- 

233 bouring Riemann problem - such as when multiple Riemann problems are solved on a 

234 sequence of grid cells in a Godunov-type numerical scheme - the local Riemann problems 

235 cannot be taken in isolation. 

236 

237 Zero rarefaction with Cr > pm and Ci — Pm- 

238 Next we consider the mollified Riemann problem with Ci = pi = pm and Cr = Pr > 

239 Pm- The convex hull and solution are depicted in Fig.jS] and the solution again consists 

240 of a contact line and zero rarefaction wave. The main difference from the previous case 

241 is that the contact line travels to the left with speed —7 instead of to the right. 

242 Note that since the location of the infiection point of fe{Pt) approaches pm as e — )■ 0, 

243 an additional zero shock of strength O(e^) and speed 0(l/e) is generated. For the sake 

244 of clarity, we have not included this wave in Fig. |8] since it is of higher order than the 

245 0{e) zero rarefaction and so has negligible impact on the solution in the limit as e — 0. 

246 ^-2- Interaction Between 1- Shock and Zero Rarefaction Wave 

247 When solving the double Riemann problem, we need to determine how the zero rar- 

248 efaction wave produced at X2 interacts with the 1-wave, which we will see shortly must 

249 be a shock. The details of the interaction can be studied using the method of charac- 

250 teristics for the four cases shown in Fig. [6] Since the zero wave has speed 0(l/e), we 
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Convex Hull Solution Profile 




Figure 7: The convex hull (left, shaded) and solution profile (right) for the generation of the zero 
rarefaction wave when pi = pm and pr < Pm- 



Convex Hull 



Solution Profile 




/.'(ft.. 



i/t = /,'(ft) 



Figure 8: The convex hull (left, shaded) and solution profile (right) for the generation of the zero 
rarefaction wave when pi = pm and pr > Pm- 
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must examine the shock-zero rarefaction interaction on two different time scales of length 
0(e) and 0(1). Since Ai = 0(1) as e — in each case, the problem can be simplified 
significantly by neglecting the shock dynamics on the 0(e) time scale. We will then show 
that the 1-wave approaches a constant speed as e ^ when t = 0(1). 



Case 1: Cr,Ci < pm- 

We first consider solutions to the mollified conservation law ^ having initial con- 
ditions ( 20 ) that satisfy Cr,Ci < pm as shown in Fig. The solution in this case 
consists of three waves: a shock, a zero rarefaction wave, and a contact line. The initial 
discontinuity at X2 generates a contact line (the 2- wave) that travels at speed X2 — 1, 
along with a zero rarefaction wave. A 1-shock originates from location xi and travels 
along the trajectory x = S{t), where S(t) has yet to be determined. Suppose that the 
zero rarefaction wave intersects the 1-wave at time t*; then the speed of the 1-wave 
satisfies the Rankine-Hugoniot condition 



dsit) uci 

— -Ai(<)- 



feiPm) 



Ci 



for tf^t*c:i 0{e) 



(21) 



The time evolution of the solution for t < t* \s illustrated in Fig. £ Jb) and the corre- 
sponding plot of characteristics in the a;,i-plane is shown in Fig.[TO The characteristics 
that intersect with the 1-wave for t < t* have speed equal to 1 to the left of the 1-wave, 
and speed 0(l/e) to the right. 

At time t* , the zero rarefaction wave starts to interact with the shock x = S{t) which 
decreases the density to the right of the 1-wave as illustrated in Fig. i [c) As a result. 



the shock wave attenuates leading to an increase in the shock speed Ai. Since the zero 
wave contains pe values lying within the interval [p„i — e,pm], we can determine the 
portion of the shock-zero wave interaction that occurs on the e time scale by finding the 
range of p^ for which f'^{pe) — 0(l/e). Equation ([9| implies that 



ap,) - O I -Ty 



Pe - Pn 



275 SO that we only need to determine the range of p^ where 



r] 



Pn 



= 0(1). 



Using the formula for the moUifier ([8|), it is easy to verify that (|23|) holds when 

and 



lim (at - 1 7^ 



a. e (-1,1). 



(22) 



(23) 



(24) 



277 It is only as p^ approaches the boundary of the e-region and lim£_5.o (a^ — l) =0, that 

278 f!:(pe) — o(l/e). Therefore, the values of within the zero rarefaction wave that satisfy 

279 (24) will interact with the 1-shock on the e time scale. 

280 Because Xi(t) = 0(1), we can ignore the influence of the shock over the e time scale 

281 when e 0. Therefore, we only need to determine the interaction between the shock 

282 and the zero rarefaction on the 0(1) time scale. By finding the range of pe within the 

283 zero wave where fi(pe) = 0(1), we can show that the shock speed Ai approaches a 
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Ai = (21) 
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Figure 9; Time evolution of the double Riemann solution in Case 1 when Cr, C; < pm- A similar picture 
holds in Case 4 when Cr,Ci > pm, except that the wave directions are reversed. 
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Case 1: Cr,Ci < p„ 



Case 2a: Cr > pm and 7/(7 + 1) ^ C; < p„ 
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Case 2b: Cr > pm and C; < 7/(7 + 1) < Pn 
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Case 3: C,. < Pm. < Ci 
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0(t) 
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Figure 10: Characteristic lines and elementary waves in the x,t-plane for various cases for the double 
Riemann solution. The thick lines (black) represent shocks while thin lines (blue and green) represent 
characteristics. 
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284 constant as e — > 0. By using the relationship (22), we know that there exists a (f> such 

285 that f'^ip,) = 0(1) for aU 

= Pm - e + 0(e), (25) 
where 0(e) — o(e) and <j>{e) > 0. Next, by bounding the integral 



Veis - P,n) ds [prn - £ + 0(e) - (Prn ~ e)]ri^{p„i -€ + 0(e) - p„i), 



we know from Eq. that 



0(e)7/e(-e + 0(e)), 
O(0(e)), 

/,(pe)-p. + O(0(e)), 



for the range of pc defined by Eqs. (25). Therefore, as e — )■ 0, we have that feipe) 



287 gfiPm) and the shock speed Ai — >■ 1. This results in a solution of the form 

Ci, if a; < xi + Xit, 
p{x,t) ^ { Pm, ii xi + Xit X X2 + X2t, 
Cr, if a; > a;2 + ^2^, 



(26) 



where Ai = A2 = 1. Therefore, at longer times the solution takes the form of a "square 
wave" propagating to the right at constant speed as pictured in Fig. 

Case 2: Ci < Pm < Cr- 

Next, we consider the double Riemann problem when Cr > pm and Ci < pm which 
generates a zero rarefaction wave and contact line at a;2, and a shock at Xi. Before the 
1-shock and zero rarefaction wave interact at time t* = 0(e), the shock speed Ai satisfies 



295 the Rankine-Hugoniot condition (21 ). For t > t* , the shock and zero rarefaction interact. 



thereby causing the value of pe to the right of the 1-wave to increase. Using a similar 
argument as in Case 1, we can deduce that when t — 0(1), 



Ai = 



=(Pm) - f{Ci) 



Pr. 



-c, 



and 



A2 



(27) 



as e 0. Note that these wave speeds are consistent with the Riemann problem in 



Eq. (17 1. The time evolution of the solution is illustrated in Fig. 



11 



There are actually two distinct sub-cases that need to be considered here, correspond- 
ing to whether C; ^ 7/(7+1) (which we call Case 2a) or C; < 7/(7 + 1) (Case 2b). 
In Case 2a, the two elementary waves (1-shock and 2-contact) do not interact, while in 
Case 2b we have Ai > A2 and so the elementary waves collide to form a single shock that 



has speed given by Eq. (19). This distinction is illustrated in the characteristic plots in 
Fig. [To] 

Case 3: Cr < pm < Ci. 

We next consider the situation when Cr < Pm < Ci, where a zero rarefaction wave 
and contact line are produced at a;2 and a shock wave is generated at a;i. Before the 
1-shock and zero rarefaction wave interact at time t* — 0(e), the shock speed Ai satisfies 
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Figure 11: Time evolution of the double Riemann solution in Case 2a when C'r > Pm and 7/(7 + 1) ^ 
C; < Pm- 
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Figure 12: Time evolution of the double Riemann solution in Case 3 when Cr < Pm < C;. 
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the Rankine-Hugoniot condition (21 ) as illustrated in Figs. a) and |(b)[ Once the zero 
wave collides with the shock, the shock speed Ai increases according to the Rankine- 
Hugoniot condition since the value of to the right of the 1-wave decreases (shown in 
Fig. l^|c) ). By ignoring interactions on the 0(e) time scale, we find that 



Ai 



9f{Prn) - f{Cl) 
Pm — Cl 



and At = 1. 



(28) 



315 
316 
317 
318 



for t = 0{1) as e — > 0. Notice that the wave speeds Ai and A2 are identical to those for 
the Riemann problem considered in Case C from Section [3j 



Case 4: Cr,Ci > Pm- 

319 The final remaining case corresponds to Cr,Ci > pm and leads to a solution satisfying 
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Eq. (261 with 



Ai 



9c{p.^) - f{Cl) 



and 



A, 



-7- 



(29) 



Pm - Cl 

The speed of the 1-shock is determined simply by the Rankine-Hugoniot condition where 
f{pm) is evaluated on the congested flow branch. This should be compared to the 1-wave 
in Case 3 where the shock speed is determined by Eq. (28 1 in which f{pm) is evaluated 
on the free flow branch. No plots of the solution have been provided for Case 4 since 
the structure is essentially the same as in Case 1 except that all waves propagate in the 
opposite direction. 



327 5. Finite Volume Scheme of Godunov Type 

328 In this section, we construct a finite volume scheme of Godunov type for our discon- 

329 tinuous flux problem based on ideas originally developed by Godunov [18] for solving the 

330 Euler equations of gas dynamics. In Godunov's method, the spatial domain is divided 

331 into cells [xj-i/2j Xj+1/2] of constant width, Ax ^ 2:^+1/2 ~ a^j-i/2j and the solution is 

332 assumed to be piecewise constant on each grid cell. Cell-averaged solution values 

333 (located at cell centers) are updated using the exact solution from local Riemann prob- 

334 lems evaluated at interfaces between adjacent cells. In each time step, the following 

335 three-stage algorithm, referred to as the Reconstruct-Evolve-Average or REA algorithm 

336 in [28', is used to update the Q": 

337 • Reconstruct a piecewise constant function p{x, tn) = Q" for all x G [xj+1/2, 2^^-1/2] 

338 from the cell average Q" at time i„. 

339 • Evolve the conservation law exactly using initial data p{x,tn), thereby obtaining 

340 p{x,tn+i) at time tn+i = t„ + At. 

• Average the solution p(a;,i„+i) to obtain new cell average values 

= — / p(x,t„+i)dx. 

341 For problems with smooth flux, the evolution step can be performed by solving a local 

342 Riemann problem at each cell interface having left state pi — and right state pr = 

343 Qj+i- As long as the time step is chosen small enough, the elementary waves produced 

344 at each interface do not interact (remembering that waves travel at finite speed in the 

345 smooth case) and hence the evolution step yields an appropriate approximation of the 

346 solution. However, as we have already shown, when the flux is discontinuous the presence 

347 of zero waves travelling at inflnite speed gives rise to long-range interactions between 

348 local Riemann problems. Consequently, we make use of solutions to the double Riemann 

349 solution derived in the previous section that incorporate the effects of zero waves. We 

350 note that the resulting algorithm has some similarities to the method of Gimse [T^. 

351 We next provide details of our implementation using LeVeque's high resolution wave 

352 propagation formulation |28j , in which the Riemann solver returns a set of wave strengths 

353 y^j-^-lf2 ^^'^ speeds s^_|.i/2 generated at each interface between states Qj and Qj+i (in 
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354 what follows, we will omit the superscript n when it is clear that time index n is assumed). 

355 The evolution step of the algorithm above can then be written as 



373 
374 
375 
376 
377 
378 
379 



At 



L^X ' ' \ •• 'I - / -I 'I - 

V 



(30) 



356 where (s)"*" = niax(s,0) and (s)^ = niin(s, 0). When the flux is smooth, the wave 

357 speed Sj_(-i/2 depends only on the states Qj and Qj+i, whereas for a discontinuous flux 

358 function this is no longer the case. When Qj, Qj+i ^ Pm, the Riemann solver yields the 

359 "standard" elementary waves whose strength and speed are given in Section [3j Because 

360 our flux is strictly non-convex, we also observe compound waves which are described 

361 within Cases 1 and 2 in Section |3] 

362 Our Riemann problem solution diverges from the standard one when either Qj = p,„ 

363 or Qj+i — Pm, in which case we construct the solution to a local double Riemann 

364 problem that requires the appropriate 1- or 2-wave given in Section [4j Note that the 

365 double Riemann solution should be viewed as two separate local Riemann problems that 

366 each produce one elementary wave: the left Riemann problem corresponds to the 1-wave 

367 and the right Riemann problem corresponds to the 2-wave. 

368 For example, if Qj — /?,„, then we construct the double Riemann problem with 

369 Qj — Pm and Qj+i — Cr, thereby obtaining the strength and wave speed corresponding 

370 to the 2-wave in Section |4j Combining together all four cases in Section [42] the speed 

371 can be written as 

'^+^/^ " I -7 , > ' ^^^^ 

372 which we note is independent of C;. 
Alternatively, if Qj+i = Pm then we construct the double Riemann problem with 

Qj = Ci, Qj+i — Pm and unknown Qj — C'r, thereby obtaining the strength and wave 
speed for the 1-wave in Section |4j In contrast with the case Qj = pm just considered, 
the 1- wave's speed depends on values of C;, pm, and Cr- Therefore, when determining 
the speed of the wave at the interface between Qj = Ci and Qj+i = Pm, we must look 
ahead to find the value Qj =^ Cr corresponding to the first value of Qj not equal to pm', 
that is, 

/ = min {fc|j + l<fcsCiV and Qk ^ Pm} ■ (32) 



380 In summary, when Q"^i = Pm, the wave speed reduces to 

5/(Pm) - fiQj) 



''^+^/^-\ 9c{Pm)-%,) 

Prn Qj 



Qi < Pn 
I > Pn 



which is visualized in Fig. 13 



382 Note that we have not yet discussed the two simple cases when Qj,Qj+i < pm and 

383 Qj, Qj+i > Pm, both of which reduce to the linear advection equation and can be trivially 

384 solved. A summary of wave strengths and speeds for all cases is presented in Table [T] 

385 A slight modification to the local Riemann solver is required to deal with the fact 

386 that algebraic operations are actually performed in floating-point arithmetic. It is highly 
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Figure 13: The four possible cases that can arise when calculating the shock speed in Eq. l|33| 



Table 1: Wave speed and strength for the 1- and 2- waves in all possible cases. 



Case 


„1 

Sj + 1/2 




„2 

■Sj + 1/2 




Qj ~ Qj+i ~ pm 














Qj = Pm and Qj / Qj+i 


Eq. (311 


Qj+l — Qj 








Qj+i = Pm and Qj / Qj+i 


Eq. (331 


Qj+l — Qj 








pm. < Qj and Pm < Qj+i 
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unlikely that the numerical value of Qj ever exactly equals Pm, and yet we find that it is 
necessary to employ the solution of the double Riemann problem when Qj is close to p„j. 
Therefore, we need to relax the requirement slightly for the zero-wave cases by replacing 
the condition Qj ~ p„i with 

\Q,~P^\^5, (34) 

where (5 is a small parameter that is typically assigned values on the order of 10^^. The 
choice of (5 is a balance between accuracy and efficiency in that taking a larger 5 value 
gives rise to significant deviations in the height of the pm plateau regions, and hence also 
errors in mass conservation. Taking 5 values much smaller than 10~^ does not improve the 
solution significantly but does require a smaller time step for stability. This modification 
influences the accuracy of the simulations and also creates an artificial upper bound on 
the maximum wave speed, not including zero waves. For example, consider the shock 
solution in Eq. (17), where as pi approaches pm the shock speed becomes unbounded. 
By enforcing (34), the parameter 5 determines how close pi can be to pm before the 
algorithm switches to the double Riemann solution. 

The time step At is chosen adaptively to enforce stability of our explicit update 
scheme, using a restriction based on the wave speeds from all local Riemann problems. 
In particular, we take 

At = CFL • min ( , , | , 

where < CFL < 1 is a constant chosen to be around 0.9 in practice. As long as the 
parameter 5 is not taken too small, this condition is sufficient to guarantee stability. 
Because the effects of the zero waves have been incorporated directly into the Riemann 
problem, they have no direct influence on the stability restriction. 

The Riemann solver described above forms the basis for the flrst-order Godunov 
scheme. We have also implemented a high resolution variant using wave limiters which 
limit the waves yVj^i/2 iii ^ manner similar to the limiting of fluxes in fiux-based finite 
volume schemes. The details of this modification are described in E51. 



409 6. Numerical Results 

410 We now apply the method described in the previous section to a number of test 

411 problems. For the high resolution scheme, we employ the superbee limiter function. In 

412 all cases, we use the discontinuous, piecewise linear fiux function ([3])-([4]) with parameters 

413 Pm = 0.5 and 7 = 0.5 that is pictured in Fig. [M] 

414 6.1. Riemann Initial Data 

415 As a first illustration of our numerical method, we use three sets of Riemann initial 

416 data: 

417 A. pi — 0.9 and p,. = 0.2, 

418 B. pi = 0.4 and pr ~ 0.9, 

419 C. pi — 0.3 and p,. — 0.98, 
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Figure 14; Plots of the discontinuous flux function used in our high resolution simulations (with pm = 0.5, 
7 = 0.5), the regularized flux for WENO calculations in Section |6.2[ and the continuous flux used in 
Section 16.31 



for which the exact solution can be determined using the methods in Sections [3] and 
|4j These initial data labeled A, B and C correspond to the three "Cases" with the 
same labels in Section [3] and pictured in Figs. [3] |4] and |5] respectively. We note that the 
numerical values for the left and right states in Case C differ slightly from the initial 
data used in Fig. [5] in order to generate larger wave speeds. 



In Fig. 15 (left), we compare the results from the first order Godunov method and 

426 the high resolution scheme with wave limiters. In each case, we perform a convergence 

427 study of error in the discrete LI norm for grid resolutions ranging from Aa; = 0.05 to 

428 0.0025. Convergence rates estimated using a linear least squares fit are summarized in 

429 Table [2] As expected, the numerical scheme converges to the exact solution for all three 

430 test problems. Godunov's method captures the correct speed for both the shock and 

431 contact discontinuities, although there is a more noticeable smearing of the contact line 

432 which is typical for this type of problem. The LI convergence rates are consistent with 

433 the order V Aa; spatial error estimate established analytically for discontinuous solutions 

434 of hyperbolic conservation laws having a smooth flux j^Hl 132] ■ The convergence rates 

435 in the L2 norm are also provided for comparison purposes and are signiflcantly smaller 

436 than the corresponding LI rates, as expected. 

437 6.2. Smooth Initial Data, With WENO Comparison 

438 For the next series of simulations, we use the smooth initial data 



(35) 






Figure 15: Left: Computed solutions for tlie Riemann initial data at t = 0.2 where S = 10 , Ax = 0.01, 
CFL = 0.95, pm = 0.5, and 7 = 0.5. Right: Convergence study in the LI norm. 
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Table 2: Convergence rates for the three Riemann problems pictured in Fig. |15| 



Godunov High Resolution 



Pi 


Pr 


LI 


L2 


LI 


L2 


0.9 


0.2 


0.643 


0.367 


1.022 


0.569 


0.4 


0.9 


0.488 


0.232 


0.832 


0.375 


0.3 


0.98 


0.754 


0.373 


1.053 


0.627 


0.1 


0.4 


0.487 


0.145 


0.700 


0.238 



where ct = 0.1 and the domain is the interval [—1,1] with periodic boundary conditions. 
This corresponds in the traffic flow context to a single platoon of cars travelling on a ring 
road, where the initial density peak has a maximum value of 1 and decreases smoothly 
to zero on either side. Fig. 16 depicts the time evolution of the solution computed using 
our high resolution Godunov scheme. Initially, the platoon begins to spread out and a 
horizontal plateau appears on the right side at a value of p = p„i. This plateau value 
corresponds to the "optimal" traffic density that is the maximum value of p for which cars 
can still propagate at the free-flow speed. The plateau lengthens as a shock propagates 
to the left into the upper half of the density profile, reducing the width of the peak. 
At the same time, the cars in the dense region spread out to the right as the plateau 
also extends in the same direction, while the left edge of the platoon remains essentially 
stationary until the dense peak is entirely gone. When the peak finally disappears (near 
time t — 0.18), the remaining platoon of cars propagates to the right with constant speed 
and unchanged shape. Note that the traffic density evolves such that the area under the 
solution curve remains approximately constant, since the total number of cars must be 
conserved. 

For comparison purposes, we have also performed simulations using the WENO 
scheme in CentPack [5], which employs a third-order CWENO reconstruction in space 
[551 and a third-order SSP Runge-Kutta time integrator jT^]. Because this algo- 
rithm requires the flux function to be continuous (although not smooth), we have used 
a regularized version of the flux that is piecewise linear and continuous, replacing the 
discontinuity by a steep line segment connecting the two linear pieces over a narrow 
interval of width It^ (instead of using the function /c in ([?]) because that would require 
an integral to be evaluated for every flux function evaluation). The regularized flux is 
shown in Fig. [Til 

From Fig. |16[ we observe that the WENO simulation requires a substantially smaller 
time step and grid spacing in order to obtain results that are comparable to our method. 
In particular, the WENO scheme requires a time step of At = O(10~^) and a spatial 
resolution of Ax — ©(lO""*) when e^, — lO^'^; this can be compared with our high 
resolution Godunov scheme for which we used a time step of At — O(10~^) when Ax = 
O(10~^) and 8 = 10~^. This performance difference is magnified further as decreases 
due to the ill-conditioning of the regularized-flux problem. 

In Fig.|17[ we magnify the region containing the density peak at time < = 0.1 to more 
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Figure 16: Time evolution of the solution with smooth initial conditions using a high resolution scheme 
{5 = 10-^, Ax = 0.005, CFL = 0.9) and a third-order WENO scheme g] [Ax = 0.0004, e^, = lO'^, 
CFL = 0.9) where p„ = 0.5 and 7 = 0.5. 



easily visualize the difference between the two sets of results for different values of 
and S. From these plots we observe that the two solutions approach one another as both 

and d are reduced, which provides further evidence that our high resolution Godunov 
scheme computes the correct solution. Note that the WENO scheme does yield a slightly 
sharper resolution of the shock than our method, but on the other hand it significantly 
underestimates the height of the plateau region when is too large. 

Next, we estimate the error in our high resolution Godunov scheme by comparing 
the computed solutions on a sequence of successively finer grids with Axp — Axq / 3^ for 
p = 0,1, . . . , P, with P = 6 levels of refinement. The finest grid solution with Ax = Axp 
is treated as the "exact" solution for the purposes of this convergence study. The errors 
a,tt — 0.05 shown in Fig. l^|^a) exhibit convergence rates of approximately 1.125 in the LI 
norm and 0.632 in the L2 norm which are consistent with the results from Section lOl We 
note that even though the initial data are smooth, our high resolution Godunov scheme 
does not obtain second order accuracy because of the shock that appears immediately 
on the right side of the plateau. 

Because convergence rates provide only a rough measure of solution error, we can 
gain additional insight into the accuracy of the method by measuring conservation error. 
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Figure 17: |(a)| WENO simulation for different Cu, for problem pictured in Fig. |16| where Ax = 0.0004 
and At = O(10~^)). |(b)| High resolution simulation for different & for problem pictured in Fig. |l6| where 
Ax = 0.005 and At = 0(10""')). 
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Figure 18: |(a)| Convergence study for the problem with smooth initial data at t = 0.05 using Axq = 0.2 
and P = 6 levels of refinement (refer to Fig. |16[ |. |(b)| Percentage of the initial vehicles lost during the 
simulation using different & values. 
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which is expressed in terms of the variation V in the total number of vehicles via 



Er = — where F " = Ax ^ Q] . 



Since we are using periodic boundary conditions and have no external sources or sinks, 
should remain approximately constant for all n. Indeed, Fig. l^[b) shows that the 



489 numerical scheme accumulates only a small conservation error over time, and that the 

490 rate of vehicles lost can be controlled by reducing S. 

491 6.3. Smooth Initial Data, With Continuous Flux Comparison 

492 In this final set of simulations, we compare solutions of the conservation law ([s]) with 

493 the discontinuous piecewise linear flux Q and the continuous piecewise linear flux 

( P, if sC p < p„, 

fc{p) = < ifp,„^p^l. (36) 



494 Both fluxes are illustrated in Fig. 14 



We begin by emphasizing that these two seemingly very different flux functions can 
still give rise to similar solutions to the Riemann problem with suitably chosen piece- 
wise constant initial data. For Cases A and C from Section |3] the convex hulls and 
the corresponding solutions are identical for the two fluxes. However, the convex hulls 
are different in Case B, where the discontinuous flux gives rise to the compound wave 
illustrated in Fig. |4] while the continuous flux generates a single shock (analogous to the 
solution in Fig. [s]). Furthermore, the continuous flux (36 1 does not give rise to any zero 



waves; instead, when either pi — p„i or p^. = pm, a single contact line is produced that 
satisfies the Rankine-Hugoniot condition. 

A clear illustration of the difference between these two fluxes is provided by comparing 
simulations on the periodic domain x G [—1,1] for smooth initial data 

p(x,0) = iexp(-^)+?, (37) 

where a = 0.1. This corresponds to the situation where there is a Gaussian-shaped 
congestion in the middle of free-flow traffic on a ring road. 



Fig. 19 depicts the time evolution of the solution computed using a high resolution 
Godunov scheme for the continuous and discontinuous fluxes. Note that the numeri- 
cal solution for the discontinuous flux requires application of the look-ahead procedure 
discussed in Section [5] while the continuous flux requires only the calculation of interac- 
tions between adjacent cells. Initially, the congested region begins to spread out and a 
horizontal plateau appears on the right side at a value oi p = pm. On the left edge of 
this plateau, a shock appears for the discontinuous flux in comparison with a gradual 
continuous variation in density for the continuous flux. 

The solutions are more drastically different when comparing the dynamics on the left 
edge of the congested region. For the discontinuous flux, a compound wave forms to the 
left of the peak, resulting in the formation of a second horizontal plateau at the maximum 
congested flow speed p = pm- Therefore, drivers on the far left will enter into a optimal 
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Figure 19: Time evolution of the solution with smooth initial conditions using a high resolution scheme 
[S = 10~^, Ax = 0.005, CFL = 0.9) for both the smooth and discontinuous flux. 



congested flow before being hit by tlie left-moving congestion wave. The continuous flux, 
on the other hand, produces only a single left-moving congestion wave. 

We conclude therefore that if the fundamental diagram is discontinuous, we expect to 
observe compound congestion waves. The leading congestion wave acts to push vehicles 
into some form of optimal congested flow, possibly a synchronized traffic state. This wave 
is followed up by a slower congestion wave. Note that we only expect these compound 
congestion waves when the upstream free flow traffic has sufficiently high density. For 
example, in Fig. 16 we observed only a single congestion shock because the traffic density 
to the left is too small to sustain this compound congestion wave. 



529 7. Conclusions 

530 In this paper, we derived analytical solutions to the Riemann problem for a hyperbolic 

531 conservation law with a piecewise linear flux function having a discontinuity at the point 

532 p — p,n- In the case when the left and right states in the Riemann initial data lie on 

533 either side of the discontinuity, the solution consists of a compound wave made up of a 

534 shock and contact line connected by a constant intermediate state at pm- In the special 

535 case when either the left or right state equals pm, the Riemann problem gives rise to a 
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536 zero rarefaction wave that propagates with infinite speed. Even though the strength of 

537 this wave is zero, it nonetheless has a significant impact on the solution structure as it 

538 interacts with other elementary waves. 

539 Our analytical results were validated using a high resolution Godunov-type scheme, 

540 based on our exact Riemann solver and implemented using the wave propagation for- 

541 malism of LeVeque i28^ . This approach builds the effect of zero waves directly into the 

542 algorithm in a way that avoids the overly stringent CFL time step constraint that might 

543 otherwise derive from the infinite speed of propagation of zero waves. We demonstrate 

544 the accuracy and efficiency of our method using several test problems, and include com- 

545 parisons with higher-order WENO simulations. 

546 We conclude with a brief discussion of several possible avenues for future research. 

547 First, a detailed convergence analysis of the algorithm would help to identify non-smooth 

548 error components that arise from the discontinuity in the flux, as well as the look-ahead 

549 procedure required to determine shock speeds in the case when pm = pi or p,.. Sec- 

550 ond, we would like to consider other nonlinear forms of the piecewise discontinuous flux 

551 function since the linearity assumption in this paper simplifies our analytical solutions 

552 considerably in that elementary waves arising from local Riemann problems consist of 

553 shocks and contact lines only. For example, it would be interesting to perform a detailed 

554 comparison with the results of Lu et al. |31) who considered a discontinuous, piecewise 

555 quadratic flux. Thirdly, we mention some preliminary computations of traffic flow using 

556 a cellular automaton model [37j that give rise to an apparent discontinuity in the funda- 

557 mental diagram. This connection between cellular automaton models (that only specify 

558 rules governing individual driver behaviour) and kinematic wave models (in which the 

559 two-capacity effect is incorporated explicitly via the flux function) merits further study. 

560 Finally, the situation where the flux is a discontinuous function of the spatial variable 

561 has been analyzed much more extensively (see the journal issue introduced by the article 

562 [6] , and references therein) . We would like to draw deeper connections between this work 

563 and the problem where the discontinuity appears in the density. 
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